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ABSTRACT 

, Aims. We want to investigate the growth of bar instability in stellar disks embedded in a suitable dark matter halo evolving in a fully consistent 
cosmological framework. 

Methods. We perform seven cosmological simulations to emphasise the role of both the disk-to-halo mass ratio and of the Toomre parameter, 
Q, on the evolution of the disk. We also compare our fully cosmological cases with corresponding isolated simulations where the same halo, is 
extracted from the cosmological scenario and evolved in physical coordinates. 

Results. A long living bar, lasting about 10 Gyr, appears in all our simulations. In particular, disks expected to be stable according to classical 
criteria, form indeed weak bars. We argue that such a result is due to the dynamical properties of our cosmological halo which is far from 
stability and isotropy, typical of the classical halos used in literature; it is dynamically active, endowed of substructures and infall. 
Conclusions. At least for mild self-gravitating disks, the study of the bar instability using isolated isotropic halos, in gravitational equilibrium, 
can lead to misleading results. Furthermore, the cosmological framework is needed for quantitatively investigating such an instability. 
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1 1. Introduction 

Several work on bar instability in stellar disks has been 
done using N-bo d y sph e rical halos (Sellwood ( 1981 ); 
! Athanassoula et alJ il987l): IPebattista & Sellwood I lilOOdb : 
■ Patsis & Athanassoula' ('200(A:' Maver & Wadslevl ( l2004 '). 
iCurir & M azzei ( 1999) were the first to emphasise the role of 
both the geometry and of the dynamical state of a live dark 
matter (DM) halo in enhancing the bar formation. Progressive 
efforts for improving models of th e halo have been made in the 
last ye ars jMazzei & Curirl ( )200lh : lAthanassoula & Misiriotis I 
( 1200 2 1). taking into account also the informations coming from 
the cosmological hierarchical clustering scenario of structures 
formation about density distribution and concentration of DM 
halos. In the meanwhile, the ever-growing computing power 
available to the community has made possible to start simula- 
tions of formation and evolution of galaxies in a fully cosmo- 
logical context. First works devoted to deep our understanding 
of disk galaxies in such a scenario (Governato et al. (2004); 
lAbadi el al. I (120031) ') have shown that is very difficult to ob- 
tain pure disk galaxies mainly because of the high angular mo- 
mentum loss of the gaseous component. Even with a careful 
choice of the hosting DM halo, the simulated galaxies appear to 
have o ver-ma ssive bulges compared to their disks. In a recent 
paper, iRobertson e t al. ( 2004) claim to have overcome most 
problems, however the bar instability has not yet been anal- 
Sent/ offprint requests to: A. Curir 



ysed in a cosmological framework. Furthermore, the high CPU 
cost of such simulations does not yet allow to explore the role 
of several parameters, most related to the phenomenological 
treatment of the star formation rate and feedbacks, on the mor- 
phologies of the generated galaxies (Mazzei & Curir (20q1|); 
lMazze'i1(l2003h ). 

In this work we present the first attempt to analyse the growth 
of bar instability in a fully consistent cosmological scenario. 
We embed a pure stellar disk inside a cosmological halo se- 
lected in a suitable slice of Universe and follow its evolution 
inside a cosmological framework. We want to explore how the 
bar instability behaves and what is the role of such a scenario. 
In particular we want to address, besides the role played by 
the disk-to-halo mass ratio, that of the dynamical state of such 
an halo as given by its substructure and infall, or more gener- 
ally by its evolution. Our model cannot be viewed as a general, 
"all-purpose" galaxy evolution model, since the gradual forma- 
tion and growth of the stellar disk is a fundamental component 
of the galaxy evolution itself. However our approach allows to 
vary parameters like the disk-to-halo mass ratio and the disk 
temperature, as given by Q parameter, to analyse the growth of 
the bar instability and its dependence on such parameters for 
the first time in a self-consistent cosmological framework. We 
analyse further the influence of the cosmological environment 
by comparing these results with those in an isolated scenario 
with the same halo. 

The plan of the paper is the following: Section 2 and 3 de- 
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scribe technical details, in particular the recipe for the initial 
disk + halo system and our framework, focusing on the cosmo- 
logical evolution and on the properties of the halo. In Section 
4 we present the whole set of our disk+halo simulations; in 
Section 5 we point out our results in the cosmological context 
and the comparison with isolate runs. In Section 6 is our dis- 
cussion and in Section 7 our conclusions. In the Appendix we 
analyse the robustness of our results, checking for particle res- 
olution and softening length effects. 

2. Numerical method 

Our galaxy model consists of a truncated exponential disk 
(Curir & Mazzei 1999), self-consistently embedded in a suit- 
able DM halo extracted from a cosmological simulation. To se- 
lect the DM halo, we perform a low-resolution (128^ particles) 
simulation of a "concordance" ACDM cosmological model: 
£l,n =0.3, Oa=0.7, o-g = 0.9, h=Q.7, where Q,„ is the total 
matter of the Universe, Oa the cosmological constant, erg the 
normalisation of the power spectrum, and h the value of the 
Hubble constant in units of 100 /i ' km s"' Mpc"'. The box 
size of our simulation is 25/i"' Mpc, which allows us an ade- 
quate cosmological tidal field and no boundary effects on our 
disk. The initial redshift is 20. We employ the pub lic paral- 
lel N-body treecode GADGET Igpringel et al. *2QQ1^- Our ini- 
tial c ondition code has been adap ted from the setup code of 
ART dKravtsov et al. I ( ll997l) : lKravtsovl ( ll999ll) : courtesy of A. 
Klypin). 

From this simulation we identify the DM halos at z=0 in the 
mass ' range 0.5- 5-10"/;"' Mq, with a standard friends-of- 
friends algorithm. We discard the halos belonging or near to 
overdense regions (see Sect. 3). Then we follow back the simu- 
lation and discard those which suffer significant mergers after a 
redshift of ~5. So we select one suitable DM halo with a mass 
M~10"/i"' Mq (at z= 0). We resamp l e it w ith the multi-mass 
technique described in iKlvpin et al. 1 1I2OOII1 . The particles of 
the DM halo, and those belonging to a sphere with a radius 4/1 ' 
Mpc, are followed to their Lagrangian position and resampled 
to an equivalent resolution of 1024 ^ particles. The total number 
of DM particles in the high resolution region is 1 2 1 65 1 2 which 
corresponds to a DM mass resolution of 1.21 ■ 10''/i 'Mq. The 
needed high frequency power is added without modifying the 
low-frequency Fourier phases of the CDM power spectrum in 
our low resolution run. The high resolution zone is surrounded 
by three shells with lower and lower resolution, the lowest one 
including all the remaining (not resampled) particles among the 
initial 128^ set. 

The size of the initial Lagrangian region is large enough to re- 
solve with high resolution not only the DM halo, but also its 
accreting sub-halos. The high-resolution DM halo is followed 
to the redshift z=0. We checked that no lower resolution par- 
ticles (intruders) are ever present at a radius lower than ~ 2 
/i 'Mpc from its centre. This corresponds to the particle with 
the minimum gravitational energy. 

' In the following, we will refer to the mass as the virial mass i.e. 
that enclos ed in a sphere with overdensity S = p/Pc-n = 178 ■ Of^'^ 
jNavarro et al. . 20001 



Our approach allows us to account for the cosmological tidal 
field acting on the DM halo and to accurately follow the evolu- 
tion of the selected halo in a self-consistent way. 
We carried out two sets of simulations embedding the galactic 
disk in the halo at the redshifts z = 2 and z - I respectively. 
The first choice corresponds to 10.24 Gyr down to z = in our 
chosen cosmology, the second one to 7.71 Gyr. 
Details of our mode l disk are presented elsewhere (e.g. 
ICurir & Mazzei 11199 9'). Here we summarise the main features 
of the disk. The spatial distribution of the star particles follows 
the exponential surface density law: pswrs - Poexp-(r/ro) 
where ro is the disk scale length, ro = 4/!"'kpc, and po is 
the surface central density. The disk is truncated at five scale 
lengths with a radius: R^isk - 20/z^'kpc. To obtain each disk 
particle's position according to the assumed density distribu- 
tion, we used the rejection method (Press et al. 1986). 

The vertical coordinate is extracted from a Gaussian distri- 
bution with a standard deviation equal to 1% of the disk radius. 
Circular velocities are assigned analytically to disk particles ac- 
counting for the global (disk+cosmological halo) potential, O. 
The radial velocity dispersion cr^ is assigned through a Toomre 
parameter Q. Q is initially constant at all disk radii and is it 
defined as Q - jH^, where k is the epicyclic frequency, and 
2 the surface density of the disk. According to the isothermal 
sheet approximation, the ratio of radial to vertical dispersion 
is fixed and constant through the disk, moreover the azimuthal 
dispersion is linked to the radia l dispersion via the epicyclic ap- 
proximation ('Hernauist " 1 993). The final velocity distributions 
are Gaussian, with the dispersions given above. 
Assigning a constant initial Q, we can classify easily our disks 
on the basis of the initial temperature. We explore two values 
of Q: 1.5, which corresponds to a warm disk, and 0.5, to a cold 
disk. The average Q va lue of stars in the M i lky W ay is esti- 
mated between 1 and 3 iBinnev & Tremainel ( Il987l) . however 
the evolution of such a parameter starting from high z is not 
known. 

Our model of galaxy is very simplified. Neither gas nor star 
formation are introduced since we aim to focus on the gravi- 
tational effect of the halo on the disk and to have hints on the 
gravitational feedback of the disk itself on the halo. Moreover, 
our technique is such that the CPU cost of one simulation, 
while large, is still much lo wer than the cost of galaxy for- 
mation simulation like that by Abadi el al. '('2003"), even if our 
force and mass resolution are comparable. Thus our work could 
give insights into self-consistent galaxy formation scenario. 
In the following we summarise the main steps of our approach: 

i) the halo is identified at redshift z = ; 

ii) its particles are tracked back to the selected redshift (i.e. 
z - I and z - 2), and the minimum of their potential well is 
calculated; 

iii) a sphere of radius Rspheie-^'^disk is extracted from the high 
resolution simulation; its bulk velocity and the position of its 
centre are recorded. Rsphere is chosen to ease the comparison of 
our results with previou s nu merical work on disk sta bihty, e.g. 
ICurir&Mazzeil ( ll999l) and lMazzei & Curir I ( l200ll) . Note that 

Rsphere >R.ir anA MZ> Mlf.^" ■ 

iv) the vector angular momentum ~f and the gravitational po- 
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tential <b are calculated, inside R sphere'^ 

v) the disk, in gravitational equilibrium with the potential Q> 
and rotating in a plane perpendicular to ~f , is generated; 

vi) we embed the disk in the high resolution simulation, at the 
chosen redshift, with its centre of mass in the minimum poten- 
tial well of the DM halo; 

vii) the bulk velocity of the halo is added to the star particles. 
The cosmological simulation is evolved then, in comoving co- 
ordinates, to the final redshift, z = 0. 

3. The DM halo 

After selecting the halo and resampling at the higher resolution 
the corresponding Lagrangian region, we run the DM simula- 
tion, to extract the halo properties in absence of any embed- 
ded stellar disk. The mass of our halo at z = 0, 1.03 ■ 10"/i ' 
Mq, corresponds to a radius, /?,,,> = 94.7/z"'kpc, which en- 
tails 84720 halo particles. The nearest DM halo ^ more mas- 
sive than IO'^/i^'Mq is ~ 1900/;"' kpc away from the centre of 
our halo; the less massive one, having mass of 4.6 ■ IO^/i^'M©, 
is ~ 215/i"' kpc away. Moreover, the behaviour of the density 
contrast, <5, is monotonically decreasing with the radius, and 6 
falls below the unity value at ~ 550/?"', 450/i"', 350/?"' phys- 
ical kpc away from the centre of our halo at z = 0,z = 1, and 
z = 2 respectively. Therefore, we conclude that the selected 
halo is living in an under-dense environment. 

The accretion history of our halo has been calculated as 
follows. Starting from redshift zi = 0, we identified DM ha- 
los using the pubHc halo finder SKID ^ (Stadel 2001) at the 
redshift zi - Z\ + dz\, corresponding to our previous simula- 
tion output. We then define as progenitor of our halo a SKID 
group at the redshift Z2, if at least a fraction / = 30% of its 
particles come from the halo at the redshift zi . We also identify 
as accreting field particles all the DM particles not belonging 
to any SKID group but belonging to the halo. We then iterate 
the procedure, using the simulation output corresponding to the 
redshift Z3 = Z2 + dzi; the progenitors are now all the groups 
which have at least a fraction / of particles coming from pro- 
genitors at zi or from the accreting field particles, and so on for 
the earlier redshifts. We check that the qualitative behaviour of 
the accretion history is not dependent on the value of / (we 
also tested /=20% and f-50%) and on the parameters used 
in SKID (we use a typical object size r - 5/i"'kpc but also 
the effect of 3/! ' kpc and 6/; ' kpc have been explored). From 
Fig. n we can note that the halo suffers its last major merger 
(i.e. a merger between two progenitors whose masses have a 
ratio which is not larger than 3) at z = 9. After z ~ 5, the 
most important contribution to its mass comes from accreting 
field particles. This contribution declines after z ~ 2 becoming 
less and less important. At z ~ 0.9, the total accreting mass is 
smaller than the mass of the larger sub-halo. Thus we conclude 
that or halo has no significant merger during the time hosting 
our stellar disk, nor immediately before. 

^ Halos have been identified using the friends of friends algorithm 
with a linking length Z = 0. 15, mean interparticle distances, with more 
than 8 particles. 

^ http://www-hpcc.astro.washington.edu/tools/skid.html 



The properties of the selected halo at three relevant redshifts 
are lis ten in TableU] Its dens i ty profile is wel l -fitt ed by a NFW 
form dNavarro el aL I (Il996l) : iNavarro et aP ^91^) at z < 2. 
The concentration, Cv,> here defined as Rvir/Rs, takes an high 
value, 18.1, confirming that this halo do es "form" at quite high 
redshift (see lWechsler et all (e.g. '2002^ for a discussion about 
the link between concentration and assembly history of the 
halo). The dimen sionless spin parame ter of the halo is defined 
^ ~ ^MVR fe^llP'^J^ al. lEoOll) where J is the angular 
momentum inside a sphere of radius R and V is the halo cuxular 
velocity, - GM/R. Its values in Table[l]are ne ar to the aver- 
age ones for our cosmological model (A - 0.035; lMaller et alH 
.200Z) 

4. Disk simulations 

We performed seven simulations of the disk+halo system as 
described here below (Sect. 4.1). By comparing results of such 
a set of 7 simulations (Sect. 5) with the DM-only, we disen- 
tangle the effect of the stellar disk on the halo evolution in the 
cosmological framework. Several simulations of the disk-i-halo 
isolated system are also run, to disentangle the effect of the cos- 
mological environment (Sect. 4.2 and 5.2). We used 56000 star 
particles to describe our disk; the (Plummer-equivalent) soften- 
ing length, the same for DM and star particles, is 0.5 /i ' kpc in 
comoving coordinates ^. We used a time-step criterion based on 
the local dynamical time (criterion "3" of the GADGET code), 
which provides 2 - 6 x lO"' time-steps from z = 2 to z = 
(except one case, simulation 5 of Table 2 which needs only 
~ 7000 time steps). The most CPU-expensive of our simula- 
tions needed ~ 5000 CPU hours to be completed on the SP4 
computer (CINECA computing center). 

4.1. Cosmological cases 

The main parameters and the initial properties of this set of 
simulations are listed in Table|2l 

A global stability crit eria for bar instability i n a disk galaxy 
is the one analysed in lEfstathiou et al.1 ( Il982l) . In such paper 
the parameters a r„ and f^^/^'^y/i (where v,„ is the maximum 
value of the disk rotational curve, r,„ the corresponding 
radius, a = rn~' and M is the disk mass) have been defined. 
lEfstathiou et al. I (Il982h stated the criterion — '■^'^ > 1.1 over 

" (aMC)'- 

the range 0.1 < a r„, < 1.3 for a disk model being stable to 
bar formation. The values of these parameters are reported in 
Table|2l 

Simulations cl, c2 and c3 in Table |2l refer to a cold disk 
(Q - 0.5). In simulation cl, at the final time (i.e. z = 0) the 
baryon fraction inside T?,,,,-, fh = MdisklMdisk+DM ~ 0.34, is 

* We note however that Cnfw is defined against R200, the radius 
enclosing a sphere with overdensity equal to 200 times the critical 
density of the Universe, and not against Ryj,. as here; therefore, in our 
cosmological model, it is always Cnfw < Cyir- At z = our halo has 
Cnfw ~ 14. 

' Note that, since the disk is modelled in physical coordinates and 
embedded in the cosmological halo at redshifts z = 2 and z = 1, its 
thickness is larger than the value of the Plummer softening we use. 



4 



Curir Mazzei & Murante: Bars in cosmology 




5 10 15 20 25 30 35 40 



Redshift 
Accretion history 



2 













most massive — \ — 
second most massive ---x - 
field ---^t^-- 
total accreting satellites □ 

y(: ^ 








































0' 

,-.-0'"' 










^^^^^"--^ 

^"""^--^ 

\^|_^^^_^ 

























1 2 3 4 5 



Fig. 1. Accretion history of our selected DM halo as a function of the redshift: top panel from redshift z = 40 to z = 0, bottom 
panel from z = 5 to z = 0; solid (red) line shows the mass of the most massive progenitor of the halo, long-dashed (green) line 
that of the second most massive progenitor, dotted (magenta) line the total mass of progenitors (the most massive excluded), and 
short-dashed (blue) line the total mass of field particles. 



44% less than its initial value, 0.53. The final baryon fraction 
of simulation c2 is ^ 0.16, compared with its initial value, 
0.28. Simulation c3 provides ft - 0.05 at z = 0, 50% less than 
its initial value. Simulations c4 and c5 provide the same final 
baryon's fractions as the corresponding simulations with the 
lower Toomre's parameter 

Simulation c6 and c7, which explore the role of the initial red- 
shift on the bar instability, provide quite the same final values 
of the baryon's fraction as the corresponding simulations c2 
and c3. Neither the Toomre parameter nor the initial redshift 
affect the evolution of this ratio which is driven by the mass 
of the stellar disk. While the baryon fraction of simulation cl 
is too high to be consistent with the cosmological value 0.166 
jRttori Il2003h . all the other simulations give baryon fractions 
in the allowed range. We however emphasise that the aim of 
the current work is not to build a realistic galaxy model, but 
to study the effect of different halo-to-disk mass ratios on the 
onset of the bar instability. We verify that the inclusion of the 



disk does not result in significant changes in the accretion 
history of the DM halo. 



4.2. Isolated cases 

We also performed several simulations of the isolated 
disk-i-halo system using the same halo as extracted from our 
cosmological simulations at z = 2 (Sect. 5.2 and Appendix). 
By comparing results of this set of simulations with the 
previous ones we aim to disentangle the effect of the large 
scale cosmological structure and of cosmological expansion 
on the system evolution. Moreover, such re sults are directly 
compar able both with our prev ious works dCurir & Mazzei I 
( 199 ^: iMazzei & Curir I j200ll) ') and with those in literature 
( lAthanassoula et alT 19871) . The initial and final values for 
these simulations are listed in Tables |2 and |3 We stress that 
our isolated halo is produced by a non dissipative collapse 
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Table 1. DM halo properties 



z 


c 


M,,, 


Ryi,. 


N,„rt 


A 


■^"60 


1'20 





18.0 


1.03 -10" 


94.7 


84720 


0.01 


0.60 


0.80 


1 


13.2 


7.7 -lO'" 


51.6 


63886 


0.02 


0.78 


0.92 


2 


8.3 


5.2 -lO'" 


30.9 


42660 


0.04 


0.92 


0.90 



I col: redshift 

II col: concentration parameter C,„- 

III col: virial mass, in 

IV col: virial radius, in physical kpc 

V col: number of DM particles within the virial radius 

VI col: spin parameter 

VII col: triaxiality parameter within a sphere of radius 60/?" 'physical kpc 

VIII col: triaxiality parameter within a sphere of radius 20/i 'physical kpc 



Table 2. Simulations: initial values. 



A' 


Q 




z 


Mdm 


Rdm 


ar,„ 




halo 


cl 


0.5 


1 


2 


0.64 


0.64 


1.9 


0.67 




c2 


0.5 


0.33 


2 


0.64 


1.94 


I. 


1.08 




c3 


0.5 


0.1 


2 


0.64 


6.4 


0.9 


1.68 




c4 


1.5 


0.33 


2 


0.64 


1.94 


I. 


1.08 




c5 


1.5 


0.1 


2 


0.64 


6.4 


0.9 


1.68 




c6 


0.5 


0.33 


I 


0.67 


2.0 


1.05 


1.05 




c7 


0.5 


0.1 


I 


0.67 


6.7 


I. 


1.6 




il 


1.5 


0.33 




0.64 


1.94 


1 


1.08 


cosm 


12 


1.5 


0.1 




0.64 


6.4 


0.9 


1.68 


cosm 


13 


1.5 


0.33 




0.64 


1.94 


1 


1.08 


cosm/frozen disk 


14 


1.5 


0.33 




0.95 


2.87 


0.85 


1.5 


NFW 


15 


1.5 


0.1 




0.95 


9.5 


1.27 


1.25 


NFW 



I col: simulation number and simulation type (c: cosmological simulations, i: isolated simulations) 

II col: Q initial value of the disk 

III col: mass of the disk in 5.9 x 10'" MO 

IV col: initial redshift (for the cosmological cases) 

V col: initial DM mass inside the disk radius 

VI col: initial hal o-to-disk mas s ratio inside the disk radius 

VII and VIII cols: lEfstathiou et al. ( 1982) parameters, where: a = tq"', v,„ is the maximum rotational velocity, and r,„ the corresponding radius. 
IX col: type of halo used (for the isolated cases) 



in a cosmological scenario. As a consequence its mass 
distribution is not spherically symmetric. Moreover such a 
halo is anisotropic and endowed with a spin parameter and 
substructure. Therefore it is different from the standard isolated 
halos used in literature to study bar instabilities, since it keeps 
a relic cosmological signature. When the halo is extracted 
from its cosmological environment, the large scale structure, 
the continuing matter infall and the expansion of the Universe 
no longer influence its evolution. Such a halo cannot be in 
gravitational equilibrium, because its evolution has not yet 
completed neither at z = 2 nor at z = 1, as shown in Fig.^ 
For this reason, the results presented below, concerning the 
behaviour of the disk embedded in such a "isolated" halo have 
to be compared with similar cases, i n which non-equ ihbrium 
DM ha los are u sed, a s e. g. in ICurir & Mazzei I ( fl999) : 
[Mazzei & Curir I ( l200ll) . After subtracting die CM velocity 
and embedding the disk, as described in Sect. 2 items i-v, we 
integrated the system in physical coordinates (the effect of the 
cosmological expansion is therefore ruled out in these models) 
. A further difference is that the softening length is now in 



physical units. We have at least 10000 time steps from the 
initial time to f = 10.24 Gyr corresponding to redshift 0. 
Finally, in order to disentangle the effect of the geometry 
and of the spin of an isolated halo we also performed two 
simulations using a Navarro Frenk and White (NFW) halo 
having the same virial radius and mass as our cosmological 
one. The initial and final values of these two simulations are 
listed in the two last lines of Tables|2and|3 



5. Results 

In this section we present the evolution of isodensity contours 
of the different cosmological simulations. From these con- 
tours we evaluate the final bar strengths which are reported 
in Table |3] Spatial resolution of the maps is always 0.5/?"' 
physical kpc and the box size is 40 times the spatial resolu- 
tion. Contours are computed at 11 fixed levels ranging from 
2 X lO"** to 0.015 in term of fraction of stars/spatial resolution 
within the total number-density of stars in the map. Following 
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Table 3. Simulations: final results 



A' 


Mdm 


Rdm 


S m 


Q< 




bulge 


bars in bars 


cl 


0.79 


0.8 


0.42 


0.38 


7 


y 


n 


cz 


U.I / 


L.iy 




U.44 


o 
5 


y 


n 


c3 


0.73 


7.41 


0.8 


0.07 


3.8 


n 


y 


c4 


0.78 


2.40 


0.48 


0.37 


5 


weak 


n 


c5 


0.73 


7.41 


0.70 


0.08 


6.5 


n 


y 


c6 


0.79 


2.43 


0.35 


0.42 


6.8 


weak 


n 


c7 


0.77 


7.73 


0.58 


0.16 


5.0 


n 


y 


il 


0.73 


2.21 


0.25 


0.4 


9.5 


y 


n 


i2 


0.51 


5.1 


0.68 


0.1 


8 


n 


y 


13 


0.7 


2.12 


0.3 


0.39 


10 


y 


n 


14 


1. 


3.03 


0.33 


0.42 


6 


y 


n 


15 


0.95 


9.5 











n 


no bar 



I col; simulation number and simulation type 

II col: DM mass inside the disk radius in 5.9 x 10'" MO 

III col: halo-to disk mass ratio inside the disk radius 

IV col: maximum bar strength at z = 0: strong bar. Mazzei & Curir] i200ll) require 5,„ < 0.6 

V col: bar strength evaluated according to Co mbes & Sanders PigsiT stronger bar corresponds to higher values of Q, 

VI col: major axis (physical kpc) corresponding to the maximum bar strength 
Vn col: morphology of the inner region of the disk 

VIII col: peculiar features inside the disk 



ICurir & Mazzei] (Il999l) . we define, as a measure of the bar 
strength, the maximum value of the axial ratio, S „, = b/a 
(Table |3ji: a strong bar corresponds to S,,, < 0.6 or to an el- 
lipticity, e = (1 - b/a), larger than 0.4. 

5.1. Morphologies of the stellar disk in the 
cosmological framework 

Fig.s 12] 15] and 0] show the evolution of isodensity contours 
of simulations cl, c2 and c3. More massive cold disks suffer 
stronger lopsided instability (m=l) from the beginning of their 
evolution which degenerates in the m=2 instability, i.e. the bar 
instability, later on. The less massive cold disks, i. e. DM dom- 
inated cases, show a weaker m=l instability. So bar instability 
develops before than in the corresponding more massive cases 
and the disk attempts to re-arrange before the end of the sim- 
ulation. Fig.s Js] and l6] compare face-on, side-on and edge-on 
isodensity contours of simulations cl and c2 at z = 0. We point 
out that our cold intermediate mass case shows peanut-shape in 
the side-on view and bulge-like contours in the edge-on view. 
Therefore, in this case, a bulge could be mis-identified due to 
the bar feature. However in the warm analogous case (FigQ 
this feature does not arise. On the other hand, our more mas- 
sive cold disk shows edge-on isodensity contours with a less 
defined inner bulge and quite ticker disk-like contours in the 
outer regions. Its side-on view corresponds to a boxy image 
without an extreme peanut feature. 

Therefore the halo-to-disk ratio has a significant influence on 
the stellar disk at z = 0, 

The Q parameter does not influence the final (z = 0) morpholo- 
gies of our less massive disks, showing always disk-like shapes 
(FigO. This suggests to regard the cold intermediate mass case 
as a peculiar one as far as the pean ut shape, is conce rned. Such 
a feature has been recognized by ICombes & Sanders . 1.1981 ) 
as caused by vertical orbital resonances. Fig. l9] shows the iso- 



density contours of simulation c4. In this simulation the higher 
value of the Q parameter stabilises the disk against the local 
Jeans instability and the bar appears later than in the corre- 
sponding cold case (simulation c2). 

Therefore wanner disks are more stable against lopsided in- 
stability than the corresponding cold cases. Inside warmer and 
less massive disks, bars in bars, namely bar features at differ- 
ent isodensity levels, nested with twisting major axes, are also 
shown. 

Isodensity contours of simulation c6 (Fig. I10> . are quite 
similar to those of the corresponding simulation c2 (Fig. |6}, 
which however starts at z - 2. Morphologies of both simu- 
lations c6 and c7 show thinner disks than simulations c4 and 
c5 given their shorter evolutionary time (^ 7.7 Gyr instead of 
K 10.24 Gyr). 

5.1 .1 . The bar strength 

A variety of quantitative parameters have been suggested to 
evaluate the strength of the bar (see Buta & Block (2001) for 
a review). Firstly we quantify the growth of the bar instability 
studying the time evolution of the ellipticity of our isodensity 
contours as a function of their major axis, a. The strength of the 
bar depends on the density contrast accounted for, and it varies 
with the distance from the centre; different choices can change 
its value but not the trend outlined in Table|3] 
Fig- E] shows that in simulation cl the strength of the bar in- 
creases with time. The length of the bar depends on the redshift 
too: it grows until z = 0.5, and then shrinks to z = 0. 
By comparing Fig. E] and Fig. O which show the elliptic- 
ity profiles of simulations c2 and c4 respectively, we point out 
that a greater Q, in intermediate mass disks, directly reflects on 
the bar strength: a stronger local gravitational instability, cor- 
responding to a lower Q value, triggers a stronger bar (Table 
O. For the less massive disks, Q is poorly influent on the bar 
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Fig. 2. Evolution of isodensity contours of simulation cl at 11 fixed levels (see Sect. 5). The size of all the frames is 2Qh ' 
physical kpc, here and in all the following figures in which isodensity contours are shown. 



strength. In these cases, the local Jeans instability has a small 
impact on the bar formation and evolution, which is dominated 
instead by the dynamics of the DM halo. 
The embedding redshift does not have a major impact on the 
bar strength. Its more important effect is the change of the bar 
length which can be connected with the larger time span of 
simulation c2 (or c3) with respect to simulation c6 (or c7). 
Therefore, the halo evolution between z = 2 and z - I does 
not seriously affect the disk instability, at least for the cold disk 
cases. 



ICombes & Sanders! ( Il98ll) have defined the bar strength at ra- 
dius R by using the parameter: Qt - where F™*' - 
[d<i>{R, ff)/d0]„,ax is the maximum amplitude of tangential force 
at radius R and < Fk{R) >= R{d<i>o/dr) is the mean axisym- 
metric radial force, at the same R, derived from the m - 
component of the gravitational potential. We evaluated the 
components of the gravitational force on a s uitable two di- 
mensional grid using the method described bv lButa & Block I 
(i2001 ). However information provided by such approach could 
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Fig. 3. Evolution of density contours of simulation c2 as described in Fig.|2 



be affecte d by spiral arm torque s and by some asymmetry in the 
bar itself ( lButa & Block l200ll) . Nevertheless, we succeeded in 
monitoring the behaviour of such a parameter for simulations 
c2, c4 and c6 (Fig.ll5>. The cold cases end up with almost the 
same value of Qt even if their evolution starts from different 
redshifts. The warmer case, instead, maintains a smaller value 
of the bar strength during all the evolution, in agreement with 
results obtained by using ellipticity parameter. Table |3 shows 
that the final values (i.e. at z = 0) of the bar strength evaluated 
with both t hese methods are cons istent. According to the clas- 
sification o f lButa& Block] ( l200 lb . we assign class 1 to our less 



massive barred galaxies if their evolution starts from z-2 (i.e. 
simulations c3 and c5), class 2 if they evolve from z - I (i.e. 
simulation c7), and class 4 to all the other ones (i.e. simulations 
cl, c2, c4 and c6). 

5.2. Comparisons with isolated cases 

In order to investigate the role of the cosmological framework 
on the bar instability, we perform isolated simulations using 
the same halo and the same disk-to-halo mass ratios as in our 
cosmological setting (Sect. 4.2). We plug Q=1.5 as stability 



Curir Mazzei & Murante: Bars in cosmology 



9 




Fig. 5. Face-on, side-on and edge-on view of isodensity contours of simulation c 1 at z = 
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Fig. 6. Face-on, side-on and edge-on view of isodensity contours of simulation c2 at z = 




Fig. 7. Face-on, side-on and edge-on isodensity contours at z = of simulation c4 in Table|2l 




Fig. 8. Face-on, side-on and edge-on isodensity contours at z = of simulation c3 (top panels) and c5 (bottom panels) in Table|2| 



parameter in the disk. Our results show that the less massive 
disks do not show important differences as far as the bar 
feature is concerned: both the bar strength is the same and the 
same bar in bar features arise. In Fig.[^we compare the halo 
radial density profiles of simulations c5 and i2. The density 
of the halo evolving in isolation becomes initially steeper, 
then it gradually flattens in the centre. In the outer regions, 
where the support of the cosmological environment is now 
lacking, the halo is slowly losing matter toward bigger scales 
and the profile is steadily steepening. On the other hand, the 
halo evolving in the cosmological environment, continues 
to accrete mass and small substructures from larger scales. 
Such accretion is still significant up to redshift z ~ 0.5 at 
least (Fig. [Q. Even if the dynamical evolution of the halo is 
different in cosmological and isolated simulations, the bar in 
the disk does form and evolves in a similar way. Thus we 
make the hypothesis that the common features of the two 
numerical experiments, namely the dynamical evolution and 



the anisotropy of the mass distribution, are the main engine for 
the bar instability. The large scale cosmological environment 
becomes a second order effect in the less massive disks. 
However the material accreting on the halo, which has been 
cut off with the halo segregation in a isolated system, plays a 
crucial role on the degree of the disk instability if the disk is 
not completely DM dominated. We conclude that the use of 
isolated halos in gravitational equilibrium for the study of the 
bar instability can give misleading results. 
Taking into account our previou s works in such i solate d 

non cosmologi cal framework jCurir & Mazzei I ( Il999l) : 

iMazzei & Curir we derive that live unrelaxed halos 

correspond to the most "realistic" approach available to 
simplify the picture. Even if the caveat outlined above cannot 
be forgotten, the dynamical state of the halo, as outlined in our 
works for the first time, plays a fundamental role in triggering 
and fuelling such a instability. 
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Fig. 9. Evolution of isodensity contours of simulation c4 as described in Fig.|2 




Fig. 10. Face on, side-on and edge-on isodensity contours of simulation c6 at z = 
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Fig. 11. Isodensity contours, in particular face-on, side-on and edge-on views of the inner bar of simulation c7 at z = 
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Fig. 12. Ellipticity as a function of the major axis, a (in physical kpc), of simulation cl at different redshifts: z = continuous 
line, z = 0.25 long-dashed line, z - 0.5 long and short dashed line, z - 0.75 dot-dashed line, z = 1 short dashed line, z - 1 .25 
dotted line. 
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Fig. 13. Ellipticity as a function of major axis, a, of simulation c2 at different redshifts; symbols are as in Fig II 21 
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Fig. 15. Evolution of bar strength after z — I evaluated using the gravitational torque (see text) for the simulations c2 (full line), 
c4 (dotted line), and c6 (dashed line) 



In order to disentangle the role of the halo's cosmological 
features like the prolate geometry and the spin on the instabil- 
ity, and to test the resolution effect, we produced an isolated 
halo with the same virial mass, radius and number of particles 
as our cosmological halo at z = 0, but with an isotropic 
NFW radial den sitv profile. The procedure is described by 
iHern auisT '('1993V We used a rejection technique to sample the 
density profile and we then assign a velocity to each particle 
following a local Maxwellian velocity dispersion. We checked 
that after 7 Gyr of evolution, the radial density profile of the 



halo is not changed, except for the "evaporation" of some 
particles dwelling in its outskirts. We embedded then a disk 
having same mass, radius and Q as in our simulation c5 and 
c3. These two simulations are labelled in Tableland Table 
|3]as i5 and i4. According to the classical theory (Sect. 4.1), 
in simulation i5 the bar instability would be inhibited. We 
successfully reproduced this result with our live NFW halo 
(Fig. I17> . Therefore the bar instability in simulation c5 is a 
genuine effect of the cosmological evolution and there is no 
evidence for a role of numerical noise. Moreover we note 
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Fig. 16. Radial density profiles of the DM halo in simulations c5 (dashed lines) and i2 (solid lines) at redshifts z = 0, 0.5, 1.0, 1.5 
from top to bottom for simulation c5 and at the equivalent evolutionary times for simulation i2. The three couples of profiles 
below have been divided by 10^, 10'*, 10^ for clarity. We also show a NFW density profile having a concentration parameter 
c ^ 23 (dotted line), obtained as a two parameters best-fit of the density profile of simulation c5 at z = 0. Length units are in 
physical kpc. 



that in simulation i5 the reaction of the DM halo to the disk 
immersion has not triggered a long-lived bar instability (see 
Fig.EJ. 



6. Discussion 

In this work we investigate the issue of the bar instability 
in stellar exponential disks embedded in a DM halo self- 
consistently evolving in a cosmological context. We aim to 
disentangle the effect of few well-defined disk parameters on 
this instability. We run also isolated simulations using the same 
cosmological halo to analyse the effect of the whole cosmolog- 
ical framework on the results. This paper is a re-visitation in 
such a cosmological scenario of the work bylc urir & Mazzei I 
( Il999l) . To compare our results with such paper we use the Rdm 
ratio in Table |2l Critical threshold values against m - 1 and 
m = 2 instabihty for thi s ratio have been calibrated also by 
lAthanassoula et all ( Il987h . They claim for a value around 0.81 
for the Rom ratio to inhibit the lopsided instability (i.e. m = 1) 
and around 2.2 to suppress the m - 2 swing amplification insta- 
bility. Even if these values are derived in a very simply frame- 
work, i.e. a isolated spheric al and an alytical D M halo, they are 
widely used in literature ( Bottema I (t2003.) : .Elmegreen el al. I 
l|2003|)), therefore we will refer to this parameter to analyse our 
initial condition, in addition to the Efstathiou et al. parameter 
introduced in Sect. 4. 1 . 

Looking at Table |2l we point out that simulations cl, c2, c4 
and c6, which develop strong final bars (TableO, are in the in- 
stability region for both these criteria. In particular in simula- 
tion c 1 , which is below the threshold of lopsided instability too 



(lAthanassoula et al]ll987l) . the signatures of such a instability 
are clearly shown in the first phases of its evolution (Fig.|5}. On 
the other hand simulations c3, c5 and c7 are stable according to 
both the criteria above. Nevertheless a weaker bar appears and 
lasts until the end of such simulations. Therefore the classical 
parameters are not good markers of the onset of the bar instabil- 
ity. In particular, when the self-gravity of the disk is negligible, 
i.e. the disk is DM dominated, the halo structure generated by 
the cosmology plays a crucial role in triggering such a instabil- 
ity. 



Our findings agree with results of lMaver & Wadslev I ll2004l) in 
the isolated framework. They found that stellar systems with 
disk-to-halo mass ratios 0. 1 become bar unstable, regardless of 
the halo concentration and the Q value, inside halos built up 
with suitable structural parameters derived from A-CDM cos- 
mology, like their circular velocity at 7?,,,,., V,,,,., the NFW den- 
sity p rofile and the spin parameter (0.06 and 0.1). We point out 
that Mayer & Wadsley ' (2004) do not take into account cos- 
mological evolution for their halos. With the same disk-to-halo 
mass ratio Ath anassoula I ( 120021) found that such a instability is 
totally inhibited inside isotropic, non rotating halos with diff'er- 
ent density profiles (eq. 1 ofjAthanassoula 2002), in agreement 
with the result of our simulation i5 (Sect. 7). This last result, to- 
gether with those of simulations performed with different num- 
ber of disk particles and with different softening length (Sect. 
7) suggests that the development of long-living bars seen in our 
simulations is a genuine physical effect and not a numerical ar- 
tifact. 

Bar instability in the DM dominated cases is strongly af- 
fected by the halo models. Moreover structural details of the 
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Fig. 17. Evolution of isodensity contours of the simulation i5 at different evolutionary times: from top to bottom, t - 5 Gyr, 
t = 7.5 Gyr and t - 10 Gyr; xy, yz and xz projections from left to right (see text (Sect. 7) for more details). 



halo, related to the cosmological framework, drive morpho- 
logical features of the stellar disk. iMaver & Wadslevl(l2004 
found indeed a central bulge after 7 Gyr which does not ap- 
pear in our corresponding case. However such a bulge shows 
up in our intermediate self-gravitating case (simulation il in 
Table Et. This feature is also e mphasi sed in the work by 
lAthanassoula & Misiriotis I ( 120021) ') and lAthanassoula I (120031) 
for a disk-to-halo mass ratio 0.2, using the same halo presented 
in^thanassoula (2002) with the higher halo concentration. 
O ur results here are in go od qualitative agreement with those 
bv lCurir & Mazzeild 19991) concerning their simulations 3, 4, 7 
and 8, with the same disk-to-halo mass ratio as in simulation 
cl of Table |2| and also with their simulations 5 and 6 which 
correspond to a disk-to-halo mass ratio 0.2. However, we re- 
mark that their simulations 3 and 4 correspond to a relaxed 
halo, whereas 5, 6, 7 and 8 to a unrelaxed halo. In particu- 
lar the Rdm initial values of their simulations 5 and 6 are re- 
spectively above and very near the 2.2 threshold value of bar 
instability, nevertheless the bar lasts until the end of their sim- 
ulations (^ 1.5 Gyr). Their simulations 1 and 2, which corre- 
spond to a relaxed dynamical state of a halo with disk-to-halo 
mass ratio 0.2, emphasise however a very different behaviour 
as far as the bar instability is concerned: the bar forms initially 
but degenerates then in a dense nucleus. Thus we argue that 
an unrelaxed dynamical state for isolated halo systems is more 
suitable to mimic a realistic "cosmological" halo, characterised 
by evolution, substructure and in-fall. This finding is important, 
since a vast majority of the work on the bar instability assumes 
a gravitationally stable halo. 



7. Conclusions 

In this work we present the first attempt to analyse the growth 
of bar instability in a fully consistent cosmological framework. 
We investigate such a issue in stellar disks embedded in a DM 
halo self-consistently evolving in a cosmological context. We 
aim to disentangle the effect of few well-defined disk parame- 
ters on this instability. We run also isolated simulations using 
the same cosmological halo to analyse the effect of the cosmo- 
logical framework. Our results show that: 

- stellar disks of different properties, i.e. mass and Q param- 
eter, embedded in the same halo and evolving in a fully 
consistent cosmological scenario, develop long Uving bars 
lasting up to redshift 0. 

- The classical criteria to account for bar instability cannot be 
validated in a cosmological framework where a bar always 
develops, due the halo evolution. 

- The strength of the bar at z = is weakly depending on 
the Q parameter, for a given disk mass. However for the 
same disk-to-halo mass ratio, colder disks show stronger 
and longer bars. Thus the less massive warm disks entail 
the weakest bars, moreover bar in bar is a common feature 
in their face-on morphology. 

- Simulations performed embedding different disks in the 
same halo, extracted a.tz-2 from the cosmological frame- 
work, show that the effects of the large scale structures are 
negligible in the less massive, DM dominated disks. 

More over by comparing results in this work with our previous 
paper dCurir & Mazzei Il999.) . we point out that live unrelaxed 



Curir Mazzei & Murante: Bars in cosmology 



16 

halos are the more suitable approach to mimic cosmological 
halos and to analyse bar instability in the less massive disks. 
The mass anisotropy and the dynamical evolution of the DM 
halo have a crucial effect in enhancing and fuelling the bar in- 
stability, also in cases where ad hoc halo mo dels provided sta- 
bility predictions, (e.g. lAthanassoula IE003 I). The large-scale 
effects, such as the continuous matter infall on the halo and the 
infall of substructures during the whole time-span of the simu- 
lation, influence the bar strength and the details of its structure. 
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Appendix A: Numerical tests 

A.1. Initialisation procedure 

While the disk is built to be in equilibrium with the grav- 
itational potential of the halo (Sect. 2), the opposite is not 
true. This means that the halo does contract when we embed 
the disk, the amount of such a contraction depending on the 
mass of the disk. This eff'ect roughly mimics the so-called 
"adiabatic contraction" (e.g Jesseit et al. ^2002) of a DM halo 
in a self-consistent galaxy formation scenario, when the gas 
firstly collapses to form a disk. Nevertheless, in our approach 
the disk appears suddenly and also contracts in reaction to the 
augmented DM concentration. This double effect could trigger 
a bar instability. Therefore we modified the code to allow for 
the presence of a frozen stellar component. Its particles do 
produce a gravitational force on the DM ones, but they are not 
subjected to any force. We built up a isolated frozen disk-nhalo 
system as described above (Sect. 4.2). Each particle of the 
disk is initially massless; such a mass increases linearly with 
the time until a final disk-to-halo mass ratio 0.3. The system 
has been integrated in physical coordinates until T = 0.51 
Gyr, which corresponds to x 1.5 tdy„ (Sect. 4.2). After this 
period, a new disk having the same mass ratio is embedded, 
in equilibrium with the gravitational potential of the halo. 
We evolved this system up to the final time, T = 10.24 
Gyr. This simulation is presented in Table |2] and Table |3] as 
simulation 13. The presence of a very strong bar is seen during 
all the evolution. In Fig. lA. II we compare the evolution of the 
maximum bar ellipticity and of the corresponding semi-major 
axis length for simulation i3 and il, in which the disk was 
added impulsively. We conclude that the reaction of the DM 
halo to the disk immersion, and the following counter-reaction 
of the disk to the increased halo concentration, are not the 
most important ingredient in producing the long-lasting bar 
we observe. 



A.2. The mass resolution 

Using a fixed number of star particles, their mass varies by 
a factor of 10 in our simulations (Table O. This can bring to 
a different numerical scatter between star and DM particles. 
Problems connected with the resolution could be the occur- 
rence of a bar instability triggered by numerical noise, which 
causes a poor sampling of the density field in the central part of 
the halo, and the numerical he ating, which also depend s on the 
number of particles involved (iLacev &Ostriker|[l985l) . could 
artificially dump the bar instability. Therefore we run isolated 
tests with different number of star particles (A^= 9400,18800, 
56000, 560000). We outline that results are slightly depending 
on the mass resolution of disk particles (Fig. lA.2t . In particular 
also the disk with 10 times more particles than our fiducial case 
keeps the bar feature, showing a more defined bar structure 
due to the higher resolution. We conclude that, even if there 
is a slight dependence of the details of the bar structure on 
the mass resolution of star particles, our main result, i.e. the 
occurrence of a long lasting bar instability even with this 
disk-to-halo mass ratio which would be classically stable, it is 
not caused by a resolution effect. 

We are currently performing new cosmological simula- 
tions (as in Sect. 4.1) with higher mass resolution: a new 
halo consisting of 282134 particles at z = 2 (compare with 
Table ^ and a stellar disk endowed with 280000 particles. 
The disk-to-halo mass ratio, defined as in Sect. 4, is 1/10. The 
halo has a spin /l=0.026, a viral mass of 7.33 ■ 10"/i"' Mq 
at z = and a "quiet" mass accretion history, similar to that 
of the DM halo used in this work. Using such halo, evolving 
in a different cosmological environment, we aim not only to 
improve the mass resolution of simulations in the cosmological 
framework but also to check how our results are depending on 
the properties of the halo. Fig. lA. 31 shows that a bar is present 
for such a very light disk also with the higher resolution. 



A.3. The softening length 

Our fiducia l softening is similar to the one suggested by 
IPower et"aD (l2003l) in their numerical convergence analysis 
of the radial density profiles of DM halos. Accounting for 
our number of particles, they suggest a minimum softening 
length: e,„,„ * Rim/ V^200, to prevent strong discreteness 
effects. From Tabled this value corresponds to e„„„ = 0.33/z"' 
kpc at z = and e„„„ - 0.15/z"' kpc at z = 2. Thus we run 
three simulations of the same isolated disk-i-halo system as in 
simulation i2 using three different softening lengths: 0.36/?"', 
0.5/?"' and 0.65/z^'kpc. In Fig lA.4l we show the isodensity 
contours at times corresponding to redshifts z = 1 and z = in 
the cosmological framework. The bar shapes and strengths are 
similar at z = 1, where comparable inner bars (as given by two 
inner levels) are shown. We note the same similarities at z = 
with the exception of simulation with lower softening length, 
which shows more round inner contours. We conclude that, 
as far as the bar instability is concerned, varying the stellar 
softening within such range has not a significant effect on the 
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results. 



A.4. The disk radius 

The last simulation we performed was devoted to explore if the 

radius of the stellar disk can affect the triggering of bar insta- 
bility. Therefore we carried out a simulation with disk param- 
eters as in simulation il, but with both disk radius and scale 
length smaller by a factor of two. Results do not differ from 
the comparison case: the system shows a definite bar instabil- 
ity developing after 0.25 Gyr and lasting until the end of the 
simulation. 
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Curir Mazzei & Murante: Bars in cosmology 




Fig. A.2. Isodensity contours at T = 7.0 Gyr of isolated simulations i2 with different number of star particles, 9400, 1 8800, 56000 
and 560000 from top to bottom; xy, yz and xz projections are from left to right. 




Fig. A.3. xy, yz and xz projections from left to right after 7 Gyr for the high resolution cosmological run with disk-to-halo mass 
ratio 0.1. See text (Sect. 7) for more details. 
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Fig.A.4. Face-on views at T=7.71Gyr and 10.24 Gyr, corresponding to z = 1 (left) and z - (right) in the cosmological 
framework, of isolated simulations i2 with stellar softening length 0.36 /i ' , 0.5 /i ' and 0.65 /i ' kpc respectively from top to 
bottom (see text for more details) 



